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Abstract 

In this paper the exact analytical solution of the motion of a rigid body with arbitrary mass 
distribution is derived in the absence of forces or torques. The resulting expressions are cast into 
a form where the dependence of the motion on initial conditions is explicit and the equations 
governing the orientation of the body involve only real numbers. Based on these results, an 
efficient method to calculate the location and orientation of the rigid body at arbitrary times 
is presented. This implementation can be used to verify the accuracy of numerical integration 
schemes for rigid bodies, to serve as a building block for event-driven discontinuous molecular 
dynamics simulations of general rigid bodies, and for constructing symplectic integrators for 
rigid body dynamics. 

1 Introduction 

The importance of the dynamics of rigid bodies in physics and engineering has been recognized 
since the early 19th century. The early work of Euler, Hermite, Poisson, Jacobi and many others 
on exactly soluble systems lead to great advances in the field of applied mathematics as well as 
in mechanics. More recently, rigid hard spheres served as the first prototypical model for fluids 
treated by computer simulation [1]. Today, rigid bodies are used to model phenomena on a variety 
of different length scales: On molecular length scales, rigid bodies are used in the modeling of the 
microscopic dynamics of molecules in condensed phases [2-4], on a mesoscopic scale, they are used 
to construct simple models of polymers and other complex systems [5], while on a macroscopic 
level they play an important role in robotics. The dynamics of rigid bodies is also of relevance 
to the computer game industry, where many improvements in simulating rigid bodies have been 
developed [6]. Finally, on an even larger scale, many astrophysical objects such as planets, satellites 
and space crafts can be regarded as rigid bodies on certain time scales [7,8]. 

According to the laws of classical mechanics, the motion of a rigid body consists of translation of 
the center of mass of the body and rotation of the orientation of the body about its center of mass. 
In the most general case in which the body is subject to forces and torques, an analytical solution 
of the dynamics is not possible and a numerical scheme is required to integrate the equations of 
motion. With the advent of powerful computers, much work has been devoted to finding stable 
and efficient integration schemes for rotating rigid bodies [9-12]. The use of numerical integrators 
is now so widespread that they are frequently used even in cases where an analytical solution of the 
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dynamics is available. This is unfortunate, since an exact solution can not only serve as a special 
case against which numerical integration schemes may be tested, but can also be directly used in 
so-called discontinuous molecular dynamics of rigid bodies, in which the bodies perform free motion 
in between interaction events. A hard sphere gas is the prototypical example of such a system. 
In these kinds of simulations, using an exact solution of the free dynamics yields an enormous 
computational benefit compared to having to integrate the free dynamics numerically (see e.g. 
Refs. 3 and 4). Furthermore, for molecular dynamics of systems with a continuous potential, an 
exact solution is also often useful in constructing so-called symplectic integration schemes, and 
typically leads to enhanced stability [12-15]. 

Perhaps the principal reason why the analytical solution of the motion of a free rigid body is 
seldom used is that the general solution is simply not well-known. Many advanced textbooks in 
mechanics avoid discussing the motion of general rigid bodies, and only consider certain symmetric 
cases [16] in which the equations of motion are particularly simple. This probably is due to the 
fact that the general solution — apparently first found by Rueb [17] in 1834 and later completed 
by Jacobi in 1849 [18] — involves special functions, called elliptic and theta functions, which are 
perhaps less familiar than other special functions. Furthermore, even when the case of a general 
rotor is discussed [8,17-20], the treatment of its motion is often incomplete, in rather abstract form 
(using complex-valued functions) and presented in a special inertial coordinate frame rather than 
a general laboratory inertial frame. 

The goal of this paper is to demonstrate that none of these issues needs to be an obstacle to 
the use of the exact solution of the equations of motion of an asymmetric rigid body in numerical 
work. It will first be shown that although the derivation of the general solution of the equations of 
motion of a rigid body requires a bit of complex analysis, the final result can be expressed without 
any reference to complex arithmetic. Secondly, the solution will be formulated in a general inertial 
frame and for general initial conditions of the rigid body. Finally, it will be shown that the special 
functions occurring in the solution can be numerically implemented in an efficient fashion. Based 
on these considerations, the numerical implementation of the (admittedly non-trivial) motion of an 
asymmetric rigid body in the absence of forces and torques is relatively straightforward. 

The paper is organized as follows: In section [2l the motion of a rigid body is reviewed starting 
in subsection 1 2 . 1 1 with a brief overview of the properties of rigid bodies. Subsequently, in subsection 
12.21 the equations of motion are given, while in 12.31 these are specialized to the case of free motion. 
The novel part of the paper starts in section [3l where a new derivation of the exact solution of the 
time dependence of the orientation of the body is given, first in general and then explicitly for the 
spherical top, the symmetric top and the asymmetric top (sections 13. H 13.21 and 13.31 respectively). 
In sectional some further numerical issues are addressed. An example is given in section [5] and the 
paper ends in section [6] with a discussion of the results and their possible applications. 

2 Review of the motion of rigid bodies 
2.1 Rigid bodies 

The shape of a rigid body is specified by the set of all material points {fj} of the body. The points 
= {xi,yi,Zi) are three-dimensional vectors with respect to a reference coordinate frame called 
the body frame, and constitute a reference orientation of the body [20]. Furthermore, a mass rrii is 
associated with each material point. The mass distribution of the body plays an important role in 
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its dynamics. 

Since any translation or rotation acting on all points f, leaves the shape of the body unchanged, 
there is an arbitrariness in the choice of body frame which can be exploited. By a suitable trans- 
lation, one may always set the center of mass to be in the origin, i.e., 



Furthermore, using a suitable rotation, the moment of inertia tensor I can be brought to its principal 
form: 



For subsequent developments, it will be assumed that these transformations have been performed. 

Because the body frame is fixed to the body, the material points r-j are independent of time. 
However since the body itself moves through physical space, the location of each point in the 
physical coordinate system, or lab frame, must be determined as a function of time to describe its 
motion. The position of mass point i in the lab frame will be denoted by = {xi{t),yi{t), Zi(t)). 
Since the body is rigid, its motion is a time-dependent orientation-and-distance-preserving trans- 
formation from the body frame to the lab frame. The most general such transformation is a 
combination of a translation and a rotation: 



Here and below, matrix- vector and matrix-matrix products such as A^{t)ri will be denoted im- 
plicitly, i.e, without a which will only be used for inner products. For notational simplicity, 
the explicit time dependence will be omitted in most expressions, i.e., ri{t) will be denoted simply 
by Tj. Exceptions are if the time argument is equal to zero (e.g. Tj(0)) or integrated over. 

In Eq. (j2.3p . the vector R denotes the position of the center of mass at time t, while the 
orthogonal matrix represents the orientation of the body with respect to the center of mass at 
that time and is often called the attitude matrix. Note that A transforms vectors from the lab to 
the body frame, while its transpose A^ transforms vectors from the body to the lab frame. 

The motions of the different material points of a rigid body are obviously closely related. Instead 
of working with the (linear) velocities of all points, one can instead use a formulation of the dynamics 
that utilizes the angular velocity vector a; = {ui, 002,^3) of the body around its center of mass. 
This angular velocity vector is defined such that its direction coincides with the rotation axis and 
its magnitude coincides with the rate at which it rotates. As a consequence of this definition, the 
velocity of any point of the rigid body satisfies the standard relation [16] 







ri{t) = R{t)+A\t)ri. 



(2.3) 




(2.4) 



where V = R is the velocity of the center of mass. Defining the antisymmetric matrix 




(2.5) 
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and using Eq. (j2.3p . Eq. (j2.4p can also be written as 



Vi = V + W{u)A^f,. (2.6) 
Taking the time derivative of Eq. (j2.3p and using Eq. (j2.6p , one sees that w and A are related via 
A^A = W{u), (2.7) 



Given the central role of rotation matrices below, it is useful to establish some notation. A 
rotation matrix U is a special orthogonal matrix that can be specified by a rotation axis h = 
(71.1,77,2,^3) and an angle ip. Here n is a unit vector, so that one may also say that any non-unit 
vector iph can be used to specify a rotation, where its norm is equal to the angle iJj and its direction 
is along the axis n. In fact, one can express this rotation matrix as ^{ipn) = exp W{^h). The 
explicit form of this rotation matrix may be found using Rodrigues' formula [16]. 



2.2 Dynamics 

It is clear from Eq. (j2.3p that the motion of a rigid body as a function of time is determined by 
the time dependence of the center of mass vector R and the attitude matrix A. According to the 
mechanics of rigid bodies, these follow from the equations of motion 

F = P T = L. (2.8) 



Here, F is the sum of all forces acting on the body, P = AIV (with M = Y^ - rrii) is the total 
momentum, r the total torque with respect to the center of mass of the body and L = lu is the 
angular momentum. One may equivalently write 

F = MV T = -f-(I^)- (2-9) 

at 

The latter equation is more conveniently written in the body frame using u = Au:, t = At and 
I = AI A"^, which yields 

T = AA''"ia> + i£l; 

= AW{u)A^iu + ii!j 

= W{u)iu + (2.10) 



where Eq. (j2.7p was used to obtain the second equality, and the third equality was obtained using 
AW(lj) At = W(Ac^) = W(cD). (2.11) 



Writing out Eq. (|2.1Up in its components gives the so-called Euler equations. Solving these equations 
yields the time dependence of the angular velocities in the body frame. To consequently find the 
attitude matrix A, one uses Eqs. (|2.7p and (|2.1ip to find 

A = -W{u)A. (2.12) 
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2.3 Force and torque- free case 



In the special case where ah external forces and torques are zero, i.e., F = and r = 0, the 
equations of motion to solve simplify to 

V = iGj = -(:jxiu A = -W(cL>)A. (2.13) 

In the absence of forces and torques, the dynamics of the system is invariant under rotations 
and under translations in time and space. As a consequence of these symmetries, the energy E, 
momentum P and angular momentum L are conserved, where E is given by 

E = Et + Er (2.14) 

with the translational and rotational energies equal to 

Et = — li'P = -M|F|2 

\ 1-1 (2.15) 

Er = -L-1 ^ L = -u - Iu = -u - iu = -(/icD^ + /2wi + ^3^1), 

respectively. 

The time dependence of the translational part of the motion follows from the equation of motion 
V = 0, which is easily solved to obtain 

y = ^(0) .2 16) 

Since V{t) = V{0), the translational energy Et is conserved in the dynamics, which, in turn, 
implies that the rotational energy Er is also conserved. 

To solve the rotational equations of motion is less trivial. One has to solve the middle equation 
of (j2.13p . which, written out in components, reads 

hoJi = 0)20)3 (/2 - -^3) 

120)2 = 0)30)1 (/3 - /i) (2.17) 
130)3 = 0)10)2 (/i - h)- 

The general solution of this set of equations will be presented below. Given this solution, the 
remaining task is to determine the attitude matrix A from Eq. (j2.12p . 

Although the general solution for A is hard to find in any textbook, it can be found in Jacobi's 
treatise on rigid body motion from 1849 [18]. Unfortunately, for mathematical elegance, the attitude 
matrix was only specified up to a rotation at a uniform speed which was not clearly identified, so that 
its numerical implementation is not altogether clear. Furthermore, Jacobi's derivation relies heavily 
on geometric arguments and Euler angles, which often pose problems in numerical applications. 
For this reason, the next section contains a novel derivation of the exact solution of A without 
reference to Euler angles, leading to an expression which is more readily implemented. 
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3 Exact solution of the attitude matrix in the absence of torques 



While discussions of how to solve the Euler equations in the body frame are common, one rarely 
sees any mention of how to go about integrating the equation for the attitude matrix. Therefore, 
we have chosen to treat this derivation in a bit more detail than one would perhaps expect from a 
computational paper. First the general form of the solution will be derived, after which the three 
different case of rigid bodies are explicitly considered: spherical bodies (section l3.ip . symmetric 
tops (section 13. 2p , and asymmetric tops (section 13. Sp . 

To obtain the general form of the attitude matrix A, we have to solve Eq. ()2.12p . Since this is 
a linear equation, its solution can be written in the form 



A = PA(0). 



(3.1) 



Here, P is a time dependent matrix which also satisfies Eq. (j2.12p . but with initial condition 
P(0) = 1. This corresponds to a case in which the body and lab frame initially coincide, i.e., in 
which the body is in the "upright" position. In this upright lab frame, the angular momentum 
vector is given by 




L{0). 



(3.2) 



It will prove to be more convenient to work in a different lab frame, called the invariant frame 
(with vectors in this frame denoted by primed quantities), in which the angular momentum vector is 
along the z-axis and is equal to L' = (0, 0, L)^. Such a frame can be found by performing a rotation 
of the original frame through a rotation matrix T'^^(O) (the reason for the peculiar notation will 
become clear below). This rotation is not unique and can be chosen such that one first rotates 
around the z-axis until the rotated angular momentum vector no longer has a y-component, and 
then subsequently rotates around the y-axis to remove the x-component of the angular momentum 
vector as well. Denoting L± = [Lf + the combined rotation is easily shown to 



T7(0) 



/km. 

' L 



\ £±(0) 

\ L 




1 





L 



( £i(0) 
Lx{0) 
1^2(0) 
Lx(0) 





V 



1^2(0) 

Lx(0) 
^i(O) 
Lx(0) 







V 



L2{0) 
L±(0) 

L 



1^2(0)^3(0) 

LLx(O) 
^i(O) 
L±(0) 
L2{0) 
L 



L 



\ 



(3.3) 



rotation to get (0) = station to get (0) = 



This rotation is defined such that if i> is a vector in the original lab frame, then v' = 'T'^{0)v 
is the corresponding vector in the invariant frame. To find the matrix corresponding to P in 
this new frame, note that since P relates vectors to the body frame through v = Fv, so that 
V = P T'^(O) v' = P' v' , one can identify 



p' = pt;(o). 



(3.4) 



^For the special case when L±{0) = one may take T'i(O) = diag(±l, 1, ±1) where the sign is chosen according 
to L-s{Q)/L = ±1, i.e., depending on whether in the original lab frame L pointed in the positive or the negative 
z-direction. 
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Since T[(0) is a constant matrix multiplying on the right, P' also satisfies Eq. (j2.12p with initial 
condition P'(0) = T[{0). 

The convenience of this choice of frame becomes clear when P' is written in terms of its columns, 



P' 



Ui U2 U3 



and one notes that 
P' 



U3 



P' L' 

L 



L 



(3.5) 



(3.6) 



the components of which are known once the Euler equations are solved. Thus, by this choice of 
frame, one is able to determine the third column of the matrix P'. 

The remaining elements of P' can all be expressed in terms of a single time-dependent angle ip. 
This is due to the orthogonality of P', which implies that the other two columns iii and U2 must 
lie in a plane orthogonal to ii^ and must also be orthogonal to each other. Denoting ei and 62 as 
two chosen orthogonal unit vectors in this plane, one can therefore write 



ill = ei cos tp — €-2 sin if) 
U2 = ei sin ^p + e2 cos ip, 

where the unit-vectors ei and 62 are chosen to beH 



£2 



X U3 



X L 

\ez X L\ 



/ L2 N 
' Li 

V 



ei = £2 X U3 



L '/ 



(3.7) 
(3.8) 



(3.9) 



(3.10) 



where = (0, 0, 1), and we have used the fact that \ez x L\ = (Lf + £2)"^^^ = L±. Other choices of 
orthogonal unit-vectors Ci and 62 change the as-of-yet undetermined time dependent angle ■0 by a 
time independent offset. The current choice has the advantage that the matrix [ei(0) 62(0) 1*3(0)] 
is identical to T;(0), so that P'(0) = T;(0), and hence 

^(0) = 0. (3.11) 

Using Eqs. (13. 7|) and ()3.8p . P' has effectively been written as a product of two rotation matrices, 



p' = t;t'2. 



where 



ei £2 W3 



/ L1L3 

LL_i_ 
L2LS 

LLx_ 

L 



-L2. 

L± 





Lz. 

L 
L . 



(3.12) 



(3.13) 



^For the special case when L± — one may take 62 = (0, 1,0)^ and ei = 62 x W3. 
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COS ip sin 0\ 

T2=(-sinV cos^ =U(-V'2). (3.14) 
1/ 

Note that T'^ in Eq. (j3.13p is a rotation matrix because ei, 62 and U3 form an orthogonal set by 
their construction in Eqs. (|3.9p and (j3.10p . The matrix P can therefore be written as 

P = T;T'2T'if(0). (3.15) 

For the special cases of the spherical top and symmetric top, it is more convenient to express P as 
a product of two rotation matrices (implicitly also found in Ref. 8): 

P = TiT2, (3.16) 

where 

Ti = T;Ti^(0) (3.17) 

T2 = T;(0)T2T'i^(0). (3.18) 

Note that Ti and T'l are determined by the solution of the Euler equations and require no further 
manipulations once the general solution of L is known. 

Turning next to the matrix T2, noting that for any rotation R and vector x, RU(a;)R^ = 
U(Ra;) and that Ti{0) z = L{0)/L, the matrix T2 can written as a rotation by an angle of — '0 
around the axis L{0)/L, i.e., using Eqs. ()3.14p and ()3.18p . 

T2 = U(-V'L(0)/L). (3.19) 

The final task to determine P consists of deriving a differential equation for the time dependent 
angle ip and solving it subject to the initial condition ip{0) = (cf. Eq. (j3.1ip ). To obtain this 
differential equation, note that from Eqs. (|2.12p . (|2.7I) and (j3.5p . it follows that 

ill = X ui, 

and hence from Eq. (j3.7p one finds 

ei cos ip — Slip sin ip — e2 sin ^p — e2ip cos tp = —Q x ei cos tp + Q x €.2 sin ip. 

Taking the inner product with 62, using that 62 • 62 = (l/2)(i|e2p/dt = and dividing by cos^H, 
yields the differential equation for the angle, ip = VI, where the time dependent frequency can be 
expressed as 

VL = -ei 62 + 62 • (a> X ci) = -ei • 62 + cj • (ei x 62) = -ei • 62 + a> • M3 
= -ei • 62 + a> • L/L. 

It follows from Eq. and (IXTnll that 

L-i{L2Li — L1L2) 



ei • 62 



LL2 



^When cos^ = 0, one can instead take the inner product with ei and divide by sinip, with the same resuh. 
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Using Eq. (j2.17p yields Li = uj^L2 — uJ2Lz and L2 = tJiLs — ws-Li, whence = L{Liubi + L2UJ2) /L\, 
or, expressed in terms of conserved quantities and lDs only, 



The angle il) is then given by 

'4)= I dt'n{t'). (3.21) 



JO 

The formal result for the matrix P, given by Eqs. (j3.13p - (j3.19p . (j3.20p and (|3.2ip . still contains 
a time integral for ip whose integrand depends on the angular velocity component 0)3. Hence the 
time dependence of this component of the angular velocity must be specified to perform the integral 
and obtain ip. The solution of the components of the angular velocity in the body frame follows 
from the Euler equations, which can be analyzed in three separate cases depending on the values 
of the principal moments of inertia. 



3.1 Spherical top 

For a spherical top, all moments of inertia are the same: Ii = I2 = I3, and the Euler equations in 
Eq. (I2.17P become extremely simple, namely, loj = 0. For the spherical top system, all components 
of the angular velocity in the body frame are therefore constant, and hence the components of the 
angular momentum in the body frame are constant as well. It therefore follows that the matrix T'^ 
in Eq. (j3.13p remains constant in time, so that 

Ti = T;Ti"^(0) = 1. (3.22) 

The frequency can also easily be determined by noting that 0)3 is constant and Ii = I2 = I3, so 
that Eq. ((3:20]) gives 



_ L(/i^f + J2(:>i) _L _ 



leading to ip = = \uj\t. Equation (j3.19p then gives, with L{0)/L = u/Vt, 

T2 = U(-^t), (3.24) 

so that one recovers the well-known result that for a spherical rotor 

P = T1T2 = U(-c^t). (3.25) 

The rotation matrix P corresponds to a rotation by an angle of — around the axis Note 
that the minus sign in the angle arises here from the fact that if the body rotates one way, the lab 
frame, as seen from the body frame, rotates in the opposite way. 
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3.2 Symmetric top 

For the case of a symmetric top, Ii = I2 but I2 / /s- In that case, one can solve the Euler equations 
(j2.17p in terms of well-known functions: 

cJi = a)i (0) cos ijjpt + 0)2 (0) sin ojpt 

Cj2 = —tDi(O) sin + (212(0) COS (3.26) 
= ^3(0). 

Here, the precession frequency ujp is given by ujp = {1 — /3//1) lD3(0). From these equations, it is 
evident that Z_l = /i[cj^ +0)2]^/'^ is conserved for a symmetric top, which allows one to rewrite Ti 
in Eq. (j3.17p as a rotation around the z-axis by an angle of —ujpt: 

cos ujpt sin ujpt 0\ 

— sm.ujpt cosoopt =\J{—U!ptz). (3.27) 
1/ 

The final rotation angle ip can be determined by noting that 

_ LjhCol + /2(2i) _L 

Ilul + Il^l (^-^^^ 

Again this is a constant but, in contrast to the spherical case, no longer equal to |a;|, so 

il^ = nt. (3.29) 

The second rotation T2 in Eq. (j3.16p is therefore given by Eq. (j3.19p with = Lt/Ii, i.e., 

T2 = U(-i(0)t//i). (3.30) 

As a result, the rotation matrix P for a symmetric top is given by 

P = Ti T2 = U(-Wptz) U(-i(0)t//i). (3.31) 

3.3 Asymmetric top 

For a general asymmetric top, all moments of inertia are distinct: Ii ^ I2 ^ I3 ^ Ii- The moments 
of inertia can be ordered in increasing order of magnitude, and we will choose to call the middle one 
always I2, while either Ji or I3 is the largest one. In the absence of forces and torques, Eq. (j2.17p is 
integrable since the energy and the norm of the angular momentum L are conserved quantities 
in the body frame. To ensure all quantities occurring in the solutions are real-valued (rather than 
complex- valued), one needs to consider the quantities Er and L^/(2l2) and make sure that [18] 

h> h> h if Er> — 

(3.32) 

h<h< h if Er < — . 

We will refer to this as the Jacobi ordering. Either situation can always be realized by choosing 
which principal axis of the body to call the first, second or third. 
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Once the Jacobi ordering is adopted, the solution of the Euler equations is given by [18-20] 



tDi = uijn cn (^LUpt + e|m) (3.33) 

UJ2 = i02m sn (^LOpt + e|m) (3.34) 

0)3 = u^rn dn (ujpt + e|m) , (3.35) 

where sn, cn and dn are Jacobi elliptic functions [21-23], and 



Cjim = sgn((I;i(0))W ^j^^j^J^^^^ (3.36) 



Cj2m = - sgn((I;i(0)) J ^j^^j^^j'j^l (3.37) 



= sgn((i3(0))W ^j^^j^^^^^ (3.38) 



= sgn(/2 - -^3) sgn(u;3(0))W . (3.39) 

{L^-2hEn){h-h) 
^ = iL^-2HEn){h-hy 

e = -UptQ = F{Lb2{G)/i02.m\m), (3.41) 
where in the last equation, F is the incomplete elliptic integral of the first kind [21], defined a^ 

F{x\m)= / (3.42) 

To give an idea of the behavior of elliptic functions, figured] shows a plot of the elliptic functions 
for two values of the elliptic parameter, m = 0.81 and m = 0.998. It is evident from the plots that 
the elliptic functions are periodic functions of their first argument, and resemble the sine, cosine and 
the constant function 1 unless the value of m is close to one. In other words, the elliptic parameter 
m determines how closely the elliptic functions resemble their trigonometric counterparts. For 
m = 0, the functions cn, sn and dn reduce to cos, sin and 1, respectively. Note that the constant 
function 1 is reminiscent of the conservation of 0)3 (cf. Eq. (j3.26p ) in the case of the symmetric top. 
Indeed, for Ji = I2, Eq. (I3.40p shows that m = 0. 

Three more numbers can be derived from the elliptic parameter m which play an important role 
in the properties of elliptic functions. These are the quarter-period K = F(l|m), the complementary 
quarter-period K' = i^(l|l— m) and the nome q = exp(—TTK'/K), which is a parameter that appears 
in various series expansions of elliptic functions. In fact, the period of the elliptic functions cn and 
sn is equal to 4i^, while that of dn is 2K. 

Given the solution of the Euler equations in Eqs. ()3.34p - (|3.35p . the matrix T'^ in Eq. (j3.13p 
is completely specified. On the other hand, to solve for the time dependence of T2 in Eq. (j3.14p . 
the integral in Eq. (j3.2ip must be performed. Unlike to previous two cases, the integrand Q of 
'0, as given by Eq. (|3.20p . is not a constant. Despite this complication, the integral can still be 



We deviate here from the somewhat more usual notation F{a\m) = F{x\m) where x — sin a. 
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(a) m = 0.81 



(b) m = 0.998 




Figure 1: Examples of the elliptic functions cn (solid line), sn (bold dashed line), dn (dotted line) 
for (a) m = 0.81 {K = 2.28..., q = 0.10...) and (b) m = 0.998 {K = 4.50..., q = 0.33...). Also 
plotted are the cosine (short dashed line) and sine (thin short dashed line) with the same period, 
for comparison. 



performed explicitly using some properties of elliptic functions as we will now show. It will require 
a bit of complex analysis to integrate O over time, so the reader may wish to skip this technical 
part and move on to the answer in Eq. ()3.54p on page [l5] (where the function ??i is the first Jacobi 
function [21,22]). 

By definition, an elliptic function is a complex function / (of a complex variable t) that is 
doubly periodic in the sense that f{t + T) = f{t + T') = f{t), and whose singularities in the finite 
complex plane consist only of poles. The doubly periodic function is completely specified by the 
values it takes inside the fundamental parallelogram spanned by r and r' in the complex plane. 
Moreover, it may be shown that an elliptic function is determined by the singularities inside that 
parallelogram up to an additive constant [22,23]. In fact, for a elliptic function f{t) having n poles 
of order one at t = tp^^c U ~ ■ ■ ^n) in the fundamental parallelogram with residues vj, one may 
write (Ref. 22, §21.5) 



m=A, + Y,r,-lo. 



'TT{t ■ 



m 



(3.43) 



where A2 is an additive constant and '!?i(n|m) is the first Jacobi theta function [21,22]. 

In order to apply Eq. ()3.43p . the singularity structure of must be examined. A little analysis 
shows that the function dn appearing in 0)3 in Eq. ()3.35p has only two poles of order one with 
opposite residues in its fundamental parallelogram, and its periods are ujpT = 2K and ujpr' = AiK' 
(Ref. 21, §16.2). Since any rational function of elliptic functions is again an elliptic function, the 
function Q. in Eq. (j3.20p is also an elliptic function with the same periods. From the form of 
in Eq. (j3.20p . one easily sees that any pole tpoie of $7 in the complex plane is due to a zero of the 
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denominator, i.e., 



^3 {tpole 



L 



(3.44) 



Note that the poles of 0)3 itself cancel in the numerator and the denominator in Eq. ()3.20p leading 
to a limit value of L//3, and do not to lead to poles in 0,. 
Using Eq. (13351) . Eq. is solved by 



f^pipoie = ± dn 



L 



h^3r. 



m 



£ + 2Kni + 4:K'n2i, 



(3.45) 



Here rii and n2 are arbitrary integers and the two it signs are independent, thus denoting four 
possibilities. Since dn changes sign when its argument is shifted over a half-period 2K'i [21], one 
lb sign in Eq. (j3.45p can be eliminated by changing the term +4K'n2i to +2K'n2i: 



uiptpole = ± dn ^ 



L 



m 



e + 2Kni + 2K'n2i. 



Using that 

dn~^(x|m) = i \K' — F{x~^\l — m)] , 
(obtained by combining §16.3.3, §16.20.3 and §17.4.46 in Ref. 21), one finds 



'^p^pole 







1 — ni^ 









e + 2Kni + 2K'n2i. 



Noting that for negative values of cDsm, one can write K' — F{I^uj^ra/ L\l — m) 
F{I^i~b^rn/ L\l — m) + 2K' , we can rewrite this as 

^ptpole = — s + 2Kni + 2K'n2i, 

where we have defined 



m 



Sgn{LJ3m)K' - 

(3.46) 
(3.47) 



Note that -K' <r]<K'. 

The periodic structure in Eq. (|3.46p can be understood as follows. The function Q depends 
on t though dn^, and dn has periods 2K and 4K'i. Note that even though dn changes sign when 
its argument is shifted over 2K'i, this sign change leaves dn^ and hence unchanged. From these 
considerations, it is evident that the actual periods of are 2K/ujp and 2K'i/L0p. Although the 
size of the fundamental parallelogram is dictated by the function under consideration, the choice of 
its origin is free. Choosing the fundamental parallelogram to be {[—K — K'i]/ujp, [—K + K'i]/ujp), 
the two poles in the parallelogram are complex conjugates tpoie and where 



Wptpoie = -e + ii]. 
The pole structure of the function Vt{t) is illustrated in Fig. [2j 



(3.48) 
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Figure 2: Periodic pole structure in the complex plane of the function ^{t) as a function of the 
variable u = ujpt. The fundamental parallelogram is indicated by a bold dashed rectangle and the 
crosses are poles. 



It is straightforward to show that the residues of Q at these poles are given by r = | for the 



pole at tpoie and 



for the pole at t*oig, so that, using Eq. (j3.43p with r = 2K/ujp, the 



integrand Q may be written as 

i d 



2dt 



log -di {ujpt + e - irj) - log -9 i(^-^{lo pt + e + irf) rnj 



(3.49) 



and further manipulated to obtain a form that can be easily integrated, 



n = A2 + 



i d 
2di 



log 1} 1 (ujpt + e -irj) 



d 



m 



m 



[log ^1 {uJpt + e - if]) 



M- — argi&i \^—{ujpt + e - irj) 



m 



(3.50) 



The constant A2 appearing in Eq. (j3.50p can be obtained using the point t = to = —e/cop, at which 
time cD2 = (cf. Eq. (^M)) and Eq. (IHT^ gives 



nito) 



L 
h 



(3.51) 



Examining the right hand side of Eq. (|3.50|) , noting that 
dt 



log ??i (^(^p* + e -ii]) rnj 



i7]) rnj 



-KUJp 






-ir]) 


rnj 


2K 


M 




-iri) 


rnj 
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and using 



■di{u\m), Eq. (|3.49p gives 




(3.52) 



Comparing Eqs. (|3.51|) and (|3.52p . we see that 



A2 



(3.53) 



Equation (j3.50p is now readily integrated to express the angle ij) in terms of the theta function as: 



and A2 is given by Eq. (I3.53p . 

With this expression for ip, the matrix P = T'^T2Ti/f(0) is now fully specified, with Ti and 
T2 given in Eqs. (j3.13p and (j3.14p . respectively. Unlike the special cases of spherical and symmetric 
rotors, no further simplifications occur in the expressions for these matrices. 

4 Numerical issues 

In this section an efficient implementation of computing the exact solutions presented in the previous 
section is outlined. Because the rotation matrices for spherical and symmetric tops in Eqs. (j3.25p 
and (|3.3ip are easily implemented using Rodrigues' formula [16], we focus on the case of an asym- 
metric rigid rotor. For this case, the attitude matrix was expressed in terms of complex theta 
functions, which hinders a straightforward and efficient numerical implementation. First it will be 
shown how all matrix elements in the attitude matrix may expressed in terms of real quantities 
and implemented efficiently. At the end of the section, an algorithm to compute the motion of an 
asymmetric rigid body will be presented. 

4.1 Implementation 

As discussed above, the matrix P can generally be written as a product of two rotation matrices 
in the form Ti T2, or three rotations in the form T'^ T2T']^(0). For the implementation of the free 
motion of the asymmetric top, we prefer the latter because it does not require Rodrigues' formula 
to be applied twice and the matrix is sparse which allows an efficient matrix multiplication. 
The extra matrix multiplication with T'^^(O) can be circumvented by using B = T'/(0)A(0). The 
matrix T'^^ is given in Eq. (j3.13p in terms of the components of the angular momentum in the body 
frame, Lj = IjOJj- The ujj are given by Eqs. (j3.33p - (j3.35p . which involve the elliptic functions. 




(3.54) 



where 




(3.55) 
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Although the occurrence of elhptic functions may seem compUcated, there are standard numer- 
ical methods to calculate the elliptic functions sn, cn and dn [21,24,25] that are very efficient and 
which makes them no more problematic to use than standard transcendental functions such as sin 
or cos. The matrix T'^ can therefore be computed numerically in a straightforward way in terms 
of the standard elliptic functions. 

The matrix is given in Eq. (13.141) and is a rotation by an angle ip around the z-axis, with ^ 
given in (j3.54p . As is evident from Eq. (j3.14p . only sin-i/; and cos'i/' rather than the angle ■0 itself 
need to be evaluated to construct Tg. Using Eq. (|3.54p and the addition formulas for cos and sin, 
these may be expressed as 



cos ip = cos(Ai + A2t) cos arg i9i + sin(Ai + A2t) sin arg t?i 
_ cos(^i + A2t) Re + sin(^i + A2t) Im -di 
~ ^(Re??i)2 + (Imi?i)2 

sin ip = sin(^i + A2t) cos arg u — cos(Ai + A2t) sin arg u = 
_ sin(y4i +A2t)Rei}i - cos(y4i + A2t)lm'di 
~ v^(Re??i)2 + (Imi?i)2 ■ 



(4.1) 



(4.2) 



where "!?! = iDi (yi^{ujpt + £ — ir/)|m) . Clearly, to compute these expressions, the real and imaginary 
parts of i?! as well as the real constants Ai and A2 must be computed. 

Noting that the function t?i has the following series expansion in the nome q (Ref. 21, §16.27.1), 



00 



i?i(n|m) = 2q^/^ ^(-l)"g"("+i) sin[(2n + l)u], (4.3) 



n=0 



the real and imaginary parts of 'di in Eq. (|4.3p for a complex argument u = -^{ujpt + s — irj), can 
be written as 



Re^i = 2gV^ V (-l)»g"("+^) cosh + ^^"^ sin + ^^"^"^^ + 

00 (4-4) 
Im^i = -2gV^ ^(-l)"g"("+^) sinh +^1)"^ eos + ^^^^.^^^ + . 

n=0 

The convergence of these series is extremely rapid due to the appearance of the In practice 

one rarely needs more than three or four terms to get to machine precision. 

Based on the series expansion of the constant Ai given in Eq. (j3.55p can be evaluated as 
follows: 

Im'!9i(2^(e — ir/)|m) 

A,=n. + arctan Re^,(^(s - .r,)H ' ^^''^ 

where n = if Re ??i > 0, n = 1 if Re i?i < and Im > 0, and n = — 1 if Re -i^i < and Im ??i < 0. 

For the constant A2, an expansion of the logarithmic derivative of "i^i can be utilized(Ref. 21, 
§16.29.1): 

■d[{u\m) ^ g2n 

„ , , — - = cottt-t-4> 7^sm2nu. (4.6) 

^ ' ' n=l 
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Using the expression for A2 in Eq. (j3.53p . where u = inr]/ {2K) is purely imaginary and noting that 
cot iu = — icothn = (e^" — l)/(e^" + 1) and siniu = isinhu = (e" — e~")/2, one can write 



.'d[{iu\m) _ e^" + 1 
'di{iu\m) e^" 



1 ^ 1 



2n 



-,2n 



(4.7) 



n=l 



Using the series expansion in Eq. ()4.7p . A2 can be evaluated from 



TTiOr, 



where 



L 

7i 2^ 



ill 



n=l 



,2n 



^2n 



(r-r") 



(4.8) 



(4.^ 



The series in Eq. (j4.8p converges if < 1. Because —K < rj < K' and q = exp{—TTK' /K), one 
has £,q^ < q < 1 so this series converges, and, because q is typically small, usually quickly. 

The above derivation assumed that the Jacobi ordering of Eq. (j3.32p was satisfied. Although 
this can always be realized by choosing which principal axis of the body to call the first, second 
or third, as is evident from Eq. (j3.32p . the choice of principal axes depends on the initial values of 
the angular velocities. Often, instead of choosing the axis depending on the initial conditions, it is 
preferable to work with a fixed convention in which the principal axes are oriented in a particular 
way with respect to the masses of the body. In that case, one can adopt the Jacobi ordering 
convention by introducing internal variables which differ when necessary from the physical ones by 
a rotation. For instance, if Ii > I2 > I3 but it should be Ii < I2 < I3 according to Eq. (I3.32p . one 
can apply the rotation matrix 



U* 




(4.10) 



which transforms between the internal and physical choices of principal axes by exchanging the x 
and z components and reversing the y component. If the order of the moments of inertia already 
follows the Jacobi ordering, U* is effectively the identity matrix. 



4.2 Algorithm for the asymmetric top 

Based on the above, the following algorithm can be set up to calculate the position of a rigid body 
at any arbitrary time, given a set of initial conditions. For efficiency, the algorithm consists of 
two steps: an initialization routine, in which some expressions are pre-calculated, and an evolution 
routine that calculates a> and A at time t. 

In the algorithm below, the functions sn, cn, dn and F are assumed to be available, but not the 
"i^i function. See Refs. 24,25 for implementations of sn, cn, dn and F. 

The initialization routine takes an initial angular velocity vector in the body frame a)(0) = 
{uJxo,ujyo,LL!zo) and inertial moments Ix, ly and 1^ (where it is assumed that ly is the middle one), 
and pre-computes a few variables as follows (in pseudo-code, in order to facilitate implementations 
in different programming languages): 
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COMPUTE L2 = J2 ^ j2 ~ 2^ ^ j2^2^ 

COMPUTE 2S = I^tD^o + /j^(I,2g + J^cD^o 

IF (2£; > L^/Ij, AND 4 < 4) OR <i2E < L^/Iy AND 4 > 4) THEN 
SET orderflag 

SET Ii=I^. I2 = ly and h = 4 

SET cjio = Wzo , = —^yo = a;a;o 

COMPUTE L_L = + -^2^^20 



ELSE 

UNSET orderflag 
SET /i = 4 , 4 = 4 and 4 = 4 
SET wio = ujxo > ^20 = c<;yo and uso = <^20 
COMPUTE L_L = \/IiUj'io + II^Iq 



END IF 

COMPUTE B = [T'^'''(0) U*] A(0) 

COMPUTE uJi^ = sgn(a;io) " 2Eh)/{h{h - 4)) 

COMPUTE a;2m = - sgn(wio) ^(^2 - 2^4)/(4(4 - 4)) 

COMPUTE wam = sgn(a;3o) V(^^ - 2Eh)/{h{h - h)) 

COMPUTE Up = sgn(4 - 4) sgn{L03o) {L"^ - 2^4) (4 - 4)/(444) 

COMPUTE m = (L2 _ 2E4)(4 - 4)/((i^^ - 2£;4)(4 - ^2)) 

COMPUTE £ = F{uJ2o/u}2m\'m) 

COMPUTE K = F{l\m) 

COMPUTE K' = F(l|l - m) 

COMPUTE q = exp{-TrK'/K) 

COMPUTE ri = sgn(LJ3o)K' - F{huJzjL, 1 - m) 

COMPUTE 4 = exp(7r77/K) 

SET ^2 = L/h + 7ra;p(^ + l)/(2i^(C - 1)) 

SET n = 1 

REPEAT 

COMPUTE 5A2 = -(7ru;p/i^)(g2"/(l - q^'')){C - T") 

INCREMENT A2 by 5A2 

INCREMENT n by 1 
UNTIL 5A2 > machine precision 
COMPUTE iVr = log (machine precision)/ log g 
SET ro = and iq = 
FOR n = TO NT DO 

COMPUTE Cr[n] = (-l)»2g"("+i)+V4 cosh (2"+^)^'? 

COMPUTE c^n] = (-l)"+i2g''(''+^)+V4 ginh (2n+i)^>7 



COMPUTE [T']^(0)U*] 




COMPUTE [T']^(0)U*] 
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INCREMENT ro by Cr[n] sin 

INCREMENT iq by q [n] cos ^^^^^IT^ 
END FOR 

IF ro > THEN k = ELSE k = sgn(io) 
COMPUTE Ai = arctan(io/ro) + kn 

STORE orderflag, h, h, h, LOim-,^2m,^-5m.,'m,ujp,e, Ai, A2, NT,Cr[\,Ci[\ and B 
END Initialization 

The evolution routine can use the pre-computed expressions in the following way: 
Evolution(t) 

COMPUTE uji = ujim cn{ujpt + e\m) 

COMPUTE UJ2 = uJ2mSn{ujpt + e|m) 

COMPUTE UI3 = uJsm dn{u;pt + e\m) 

SET Re ??i = and Im i?i = 

FOR n = TO NT DO 

INCREMENT Rei^i by Cr[n]sm{{2n + l)7r{u}pt + £)/{2K)) 
INCREMENT liwdi by Ci[n]cos{{2n + l)Tr{ujpt + £)/{2K)) 

END FOR 

COMPUTE C = cos(^i + A2t) , S = sin(yl i + A2t) 
COMPUTE cosV' = (C Rei?i + 5 Im ) / y ^Re i^f + Im 'dj 
COMPUTE sinV' = (S Re-di - C Im ) / y^Re i?f + Im 
COMPUTE L_L = y^I'iio'i + I'iu'j 
IF orderflag IS SET THEN 

COMPUTE [U*T'J 

SWAP uji and 0)3 
CHANGE SIGN of UJ2 
ELSE 

COMPUTE [U*T[ 

END IF 

COMPUTE A = [U* T'J 

RETURN uii, UI2, 
END Evolution 

Some remarks about this pseudo-code: 

• The orderflag indicates whether the Jacobi ordering convention Eq. (j3.32p is satisfied. If 
not, U* in Eq. (j4.10p is used. If Eq. (j3.32p is satisfied, U* is set equal to the identity matrix. 

• With these definitions, P is replaced by U*PU*. This is accomplished simply by T'^^ — > 

u* t;. 

. As a consequence, U* in Eq. (j4.10p is only implicitly used in the combination [U* T^] 
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Figure 3: Example of the exact solution, (a) u}{t), (b) position of the 'axes'. For details see text. 



• The initial value of A(0) occurs as A{t) = U* T'^ T'^ T'I{0) U* A(0), so in the initiaUzation 
routine, we only need to store the combination B = T'^^(O) U* A(0). 

• The machine precision depends on the floating point precision used in the calculation; for 
64 bit double precision, this is of the order of 10~^^-10~^^. 

• For clarity of the algorithm, matrix products have not been explicitly written out, and effi- 
ciency improvements such as computing intermediate expressions and using a recursive eval- 
uation for the sin's and cos's have not been shown here. 

• An implementation of this code in C, which includes these improvements, can be found 
on the internet at http://www.chem.utoronto.ca/staff/JMS/rigidrotor.html and in the 
appendix. 



5 Example 

As an example, consider an object composed of six point masses which are arranged, in the body 
frame, at the points (a, 0,0), (—a, 0,0), (0,6,0), (0,-6,0), (0,0, c) and (0, 0, — c). All points are 
assumed to have unit mass, so that the inertial moments are given by Ii = 2(6^ -|-c^), I2 = 2(a^ -|-c^) 
and I3 = 2(a^ + 6^). Choosing a > b > c ensures Ii < I2 < h- 

In particular we will consider o = 3, 6 = 2 and c = 1, yielding Ii = 10, I2 = 20 and I3 = 26. 
As initial conditions we will take A(0) = 1 and 4*^(0) = (1,15,1). Because of the large value of 
the y-component of the angular velocity, one may expect motion to consist primarily of rotation 
around the y-axis. In figure [3^ the components of the angular velocity in the body frame have 
been plotted as a function of time, while in figure [3)3 the projection on the x — z plane of the point 
masses that in the body frame are located at (a, 0,0) (the 'long axis'), (0,6,0) (the 'middle' axis') 
and (0, 0, c) (the 'short axis') are plotted. It is clear that the rotational motion does not consist 
of small perturbations to a rotation around the y-axis. We stress once more that, up to machine 
precision, the results in figure [3] are exact. 



20 



6 Discussion 



In this paper, the general solution of the rotational motion of a rigid body in the absence of 
external torques and forces was derived. Explicit expressions for the angular velocities and the 

attitude matrix were obtained in terms of real quantities to facilitate numerical evaluation. Note 
that even though the solution of rotational motion for bodies without a simple mass distribution 
contains generalizations of the familiar sine and cosine functions, the motion typically appears quite 
complex and notably different from that of a spherical top. 

The general solution of the equations governing rigid body dynamics in the absence of forces 
and torques presented here is potentially useful in several important applications. The primary 
advantage of having in hand analytical solutions of the equations of motion of a system lies in 
the fact that all relevant properties of a soluble system can be determined at arbitrarily many 
and arbitrarily distant moments in time. Applications in which knowing exactly the position and 
orientation of a body at specific moments in time is paramount may benefit from the results 
presented here. Such applications are abundant in a wide variety of contexts. For example, in 
astrophysics, many objects such as space crafts, asteroids, certain planets and moons, behave 
on short time scales as rigid bodies. These bodies are not free since they feel typically weak 
gravitational fields. However, if their dimensions are small enough compared to the gradients of 
the gravitational field, gravitational forces effectively influence motion only of the center of mass, 
while the rotational motion is that of a free rotor described here. 

Another obvious application of the solution detailed in this paper is as a diagnostic tool for nu- 
merical integration techniques designed for rotating bodies with external torques. Such techniques 
are of considerable interest, but to establish their accuracy, one needs to be able to compare results 
of approximate integration schemes with exact results. To date most comparisons are carried out 
for free systems with a high degree of symmetry and simple rotational motion [10]. Given the 
relative complexity of motion in the asymmetric case compared with that of a spherical rotor, such 
comparisons do not appear to be very stringent. 

The exact solution is also of practical use in symplectic integrators for use in continuous molec- 
ular dynamics [15]. There are already various symplectic integrators using the exact solution of 
some part of the dynamics [12-14], which generally seems to improve stability and accuracy over 
simple splitting methods [9-11,26]. However, these do not use the exact solution of the attitude 
matrix. The integrator of Celledoni and Safstrom, for example, uses the exact solution of the Euler 
equations but uses an approximate expression for the attitude matrix [12]. Using the exact solution 
of the attitude matrix further improves the stability and accuracy of this integrator [15]. 

Perhaps the most direct application of the implementation of the exact solution is in simu- 
lating complex rigid molecular systems using discontinuous molecular dynamics methods. In this 
approach, various components of the system interact via discontinuous potentials, leading to im- 
pulsive forces and torques that act on molecules at specific moments in time [1,2,5,27,28]. As 
a result, the motion of all bodies in the system is free between impulsive events that alter the 
trajectory of the body via discontinuous jumps in the momenta or angular velocities at discrete 
"collision" times. In order to determine the time at which molecules in the simulation interact, the 
exact location and orientation of all bodies in the system must be computable at arbitrary times. 
If the configurations of the system are computed through numerical integration (using e.g. one 
of the integrators in Refs. 9-12), such simulations would become inefficient. For this reason, to 
date, most simulations of rigid bodies interacting via discontinuous potentials have been restricted 
to systems in which rotational motion is governed by the equations of a spherical rotor (see Refs. 
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2). Armed with the results of this paper, the technique of discontinuous molecular dynamics can 
now be applied to any rigid model — symmetric or asymmetric — with discontinuous interactions of 
step-potential form. Examples of such studies can be found in Refs. 4. 
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Appendix: Implementation in C 

#define MACHPREC 3.E-16 

#define MY_PI_2 1.57079632679489661923 

/* The following macro performs an in-situ matrix multiplication: */ 
#define RIGHTMULTMATRIX(M1 , M2, x, y)\ 
X = M1[0] [0] ; y = Ml [0] [1] ;\ 

M1[0] [0] *= M2[0] [0] ; Ml [0] [0] += y*M2[l] [0]+Ml[0] [2]*M2[2] [0] ;\ 
M1[0][1] *= M2[l][l]; M1[0][1] += x*M2 [0] [1] +M1 [0] [2] *M2 [2] [1] ;\ 
M2[2] [2] ; M1[0] [2] += x*M2 [0] [2]+y*M2[l] [2] ;\ 
X = Ml[l] [0] ; y = Ml [1] [1] ;\ 

M1[1][0] *= M2[0][0]; Ml [1] [0] += y*M2 [1] [0] +M1 [1] [2] *M2 [2] [0] ; \ 
Ml[l][l] *= M2[l][l]; Ml[l][l] += r*M2[0] [1]+M1[1] [2]*M2[2] [1] ;\ 
Ml[l] [2] *= M2[2] [2] ; Ml [1] [2] += r*M2 [0] [2]+y*M2[l] [2] ;\ 

X = Ml [2] [0] ; y = Ml [2] [1] ;\ 

Ml [2] [0] *= M2[0] [0] ; Ml [2] [0] += y*M2[l] [0]+Ml[2] [2]*M2[2] [0] ;\ 
Ml[2][l] *= M2[l][l]; Ml[2][l] += x*M2 [0] [1] +M1 [2] [2] *M2 [2] [1] ;\ 
Ml [2] [2] *= M2[2] [2] ; Ml [2] [2] += x*M2[0] [2]+y*M2[l] [2] ; 

/* A number of precomputed expressions are taken together in the 

following structure: */ 
typedef struct { 
int orderFlag, NT; 

double II, 12, 13, omegalm, omega2m, omega3m, omegap, freq, epsilon, 
Al, A2, L, m, *cr, *ci, B[3][3]; 

} Top; 

/* The following function performs some pre-calculations . Input: Ix, * 

* ly and Iz are the three principal inertial moments (it is assumed * 

* that Ix<Iy<Iz or Ix>Iy>Iz) , omegax, omegay and omegaz are the * 

* angular velocities in the principal axis (i.e. body) frame, and A * 

* is the initial attitude matrix. Returns a Top structure containing * 

* all necessary precomputed parameters to efficiently calculate the * 

* omega's and A's at later times (see the Evolution fimction below): */ 
Top Initialization(double Ix, double ly, double Iz, 
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double omegax, double omegay, double omegaz, double A [3] [3]) 

{ 

double 

a = Ix*omegax, /* LI, angular momentum in x-direction */ 

b = Iy*omegay, /* L2, angular momentum in y-direction */ 

c = Iz*omegaz, /* L3, angular momentum in z-direction */ 

d = a*a+b*b+c*c, /* L.L, norm squared of angular momentum */ 

e = a*omegax+b*omegay+c*omegaz, /* 2 E */ 
f, omegal, omega2, omegaS, Kp, q, r, i, s, g, h; 

int n; 

Top R; /* will contain the precomputed numbers */ 

R.L = sqrt(d) ; /* \L\ */ 

if ( (e>d/Iy && Ix<Iz) | | (e<d/Iy && Ix>Iz) ) { 
/* Check if Jacobi-ordering is obeyed: */ 

R.orderFlag = 1;/* Jacobi ordering ensured by swapping the order of x and z */ 
R.Il = Iz; /* directions and reversing the y direction */ 

R.I2 = ly; 
R.I3 = Ix; 
omegal = omegaz; 
omega2 = -omegay; 
omegaS = omegax; 
f = hypotCb, c) ; 

/* Fill the matrix R.B with transpose of T1(0) , using ordering: */ 
R.B[0][0] = -f/R.L; R.B[0][1] = b*a/R.L/f; R.B[0][2] = a*c/R.L/f; 
R.B[1][0] = 0; R.B[1][1] = -c/f; R.B[1][2] = b/f ; 

R.B[2][0] = a/R.L; R.B[2][1] = b/R.L; R.B[2][2] = c/R.L; 

} else { 

R.orderFlag =0; /* Jacobi ordering correct */ 

R.Il = Ix; 

R.I2 = ly; 

R.I3 = Iz; 

omegal = omegax; 

omega2 = omegay; 

omegaS = omegaz; 

f = hypot (a, b) ; 

/* Fill the matrix R.B with the transpose of T1(0) : */ 
R.B[0][0] = a*c/R.L/f; R.B[0][1] = b*c/R.L/f; R.B[0][2] = -f/R.L; 
R.B[1][0] = -b/f; R.B[1][1] = a/f ; R.B[1][2] = 0; 

R.B[2][0] = a/R.L; R.B[2] [1] = b/R.L; R.B[2][2] = c/R.L; 

} 

RIGHTMULTMATRIX(R.B, A, r, i) ; /* calculate B */ 

a = d-e*R.I3; /* compute four subexpressions */ 

b = d-e*R.Il; 
c = R.I1-R.I3; 
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d = R.I2-R.I3; 

R.omegalm = copysign(sqrt (a/R. Il/c) , omegal) ; 
R.oniega2iii = -copysignCsqrt (a/R. I2/d) , omegal); 
R.omegaSm = copysign(sqrt (-b/R. I3/c) , omegaS) ; 

R.omegap = d*copysign(sqrt(b/(-d)/R.Il/R.I2/R.I3) , omegaS) ; /* prec. freq */ 
R.m = a*(R.I2-R.Il)/(b*d) ; /* ellipic parameter m */ 

R.epsilon = F(omega2/R.omega2m, R.m); /* initial phase */ 

R.freq = MY_PI_2/F(1 .0, R.m); /* frequency relative to precession */ 

Kp = F(1.0, 1.0-R.m); /* K' , complementary quarter period */ 

q = exp(-2.0*R.freq*Kp) ; /* elliptic nome */ 

e = exp(R. f req* (copysign(Kp, omegaS) -F(R.I3*R.omega3m/R.L, 1.0-R.m))); 



/* = exp 



7r(sgn(a;3)j<:'-F(J3a)3^/L|m) 
2K 



*/ 



f = e*e; /* this f is equal to xi */ 

R.A2 = R.L/R.Il+R.freq*R.omegap*(f+l)/(f-l) ; /* first term in A2 series */ 

a = 1.0; /* a will be the 2n-th power of q */ 

b = 1.0; /* b will be the nth power of xi */ 

n = 1; 
do { 

a *= q*q; /* update a and b recursively */ 

b *= f; 

c = -2.0*R.freq*R.omegap*a/(l-a)*(b-l/b) ; /* the next term in A2 */ 

R.A2 += c; /* add a term of the series */ 

n++; 

} while (f abs(c/R.A2)>MACHPREC && n<10000);/* stop if converged or n too big */ 
/* determine upper bound on number of terms needed in theta function series: */ 
R.NT = (int) (log(MACHPREC)/log(q)+0.5) ; 

R.cr = (double *)malloc(sizeof (double) *(R.NT+1)) ; /* allocate memory */ 

R.ci = (double *)malloc(sizeof (double) *(R.NT+1)) ; 

a = 1.0; /* a will be the 2n-th power of q */ 

b = 1.0; /* b will be (-l)"n q-{n(n+l)} */ 

R.cr[0] = (e+l/e) ; /* zeroth term in the series for real and imag part */ 

R.ciEO] = -(e-l/e) ; 

s = sin(R. f req*R. epsilon) ; /* s = sin( (2n+l)x) , c = cos((2n+l)x) */ 

c = cos (R.freq*R. epsilon) ; 

g = 2.0*c*s; /* sin(2x) */ 

h = 2.0*c*c-1.0; /* cos(2x) */ 

r = R.cr[0]*s; /* real part of the theta function, zeroth term */ 

i = R.ci[0]*c; /* imaginary part of the theta function, zeroth term */ 

for (n = 1; n <= R.NT; n++) { 
e *= f; /* e = xi''{n+l/2} */ 

a *= q*q; /* update a and b recursively */ 

b *= -a; 

R.cr[n] = b*(e+l/e); /* compute next coefficient of theta series */ 

R.ci[n] = -b* (e-l/e) ; 
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d = s; /* compute sin and cos recursively */ 

s = h*s+g*c; 
c = h*c-g*d; 

r += R.cr[n]*s; /* add next terms */ 

i += R.ci[n]*c; 

if ( (fabs(R.cr[n]) < MACHPREC) && (f abs (R. ci [n] ) < MACHPREC) ) 
R.NT = n-1; /* if converged, adjust NT */ 

} 

R.Al = atan2(i, r) ; /* compute arg(r, i) */ 

return R; /* done, return all precomputed values */ 

} 



/* The following function calculates the angular velocities and attitude * 

* matrix at a time t . Input : R is a pointer to a Top structure containing * 

* all necessary precomputed parameters to efficiently calculate the omega's * 

* and A's, and t is a time. R should be generated by Initialize function. * 

* Output: *omegax, *omegay and *omegaz are filled with the angular * 

* velocities in the principal axis (i.e. body) frame at time t and A is * 

* filled with the attitude matrix at time t. */ 
void Evolution(Top *R, double t, double *omegax, double *omegay, double *omegaz, 

double A [3] [3]) 

{ 

int n; 

double omegal , omega2 , omegaS , r, i, u, s, c, g, h, f; 



u = R->omegap*t+R->epsilon; 



/* compute argument of elliptic functions */ 



/* compute sn, cn, dn */ 
/* multiply by the amplitudes */ 
/* of the respective ang. vel. */ 



/* check for Jacobi ordering */ 
/* if adjusted, invert x and z and reverse y */ 



SNCNDN(u, R->m, &omega2, feomegal, &omega3) ; 
omegal *= R->omegalm; 
omega2 *= R->omega2m; 
omega3 *= R->omega3m; 
if (R->orderFlag == 1) { 

*omegax = omegaS; 

*omegay = -omega2; 

*omegaz = omegal; 

omegal *= R->I1; 

omega2 *= R->I2; 

omegaS *= R->I3; 

f = hypot (omegal , omega2) ; 

/* compute Tl(t) , taking into account the ordering: */ 

A[0][0] = -f/R->L; A[0][1] = 0; A[0] [2] = omega3/R->L; 

A[l] [0] = -omega2*omega3/R->L/f ; A[l][l] = -omegal/f; A[l] [2] = -omega2/R->L; 
A[2] [0] = omegal*omega3/R->L/f ; A[2] [1] = -omega2/f; A [2] [2] = omegal/R->L; 
} else { 

*omegax = omegal; /* no adjustment necessary */ 
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*omegay = omega2; 

*oniegaz = omegaS; 

omegal *= R->I1; 

omega2 *= R->I2; 

omegaS *= R->I3; 

f = hypot (omegal , oinega2) ; 

/* compute Tl(t) : */ 

A[0][0] = omegal*omega3/R->L/f ; A[0] [1] = -omega2/f; A[0] [2] = omegal/R->L; 
A[1][0] = omega2*omega3/R->L/f ; A[l] [1] = omegal/f; A[l] [2] = omega2/R->L; 
A[2][0] = -f/R->L; A[2] [1] = 0; A [2] [2] = omega3/R->L; 

} 

s = sin(R->f req*u) ; /* s = sin( (2n+l)x) , c = cos((2n+l)x) with x = freq*u */ 
c = cos(R->f req*u) ; 

g = 2.0*c*s; /* g = sin 2x, h = cos 2x, used in recursion */ 

h = 2.0*c*c-1.0; 

r = R->cr[0]*s; /* zeroth term in series for the theta function */ 

i = R->ci[0]*c; 

for (n = 1; n <= R->NT; n++) { /* compute series */ 

u = s; 

s = h*s+g*c; /* computes sin((2n+l)x) and cos((2n+l)x recursively */ 

c = h*c-g*u; 

r += R->cr[n]*s; /* next term in series for theta function */ 

i += R->ci[n]*c; 

} 

s = sin(R->Al+R->A2*t) ; /* s = sin(A_l+A_2 t) */ 

c = cos(R->Al+R->A2*t) ; /* c = cos(A_l+A_2 t) */ 

u = s; /* use addition formula to compute s = sin psi and c = cos psi */ 

s = s*r-c*i; /* where psi = A_l+A_2 t+arg(r, i) */ 

c = c*r+u*i; 
u = hypot (r, i) ; 
s /= u; 
c /= u; 

for (n = 0; n<3; n++) { /* perform multiplication with T2' */ 

u = A[n] [0] ; 

A[n] [0] = A[n] [0]*c-A[n] [l]*s; 
A[n] [1] = A[n] [l]*c+u*s; 

} 

RIGHTMULTMATRIX(A, R->B, r, i) ; /* gives the final attitude matrix */ 

} 

Let us briefly mention a few accessible implementations of elliptic functions and elliptic integrals: 
• Using the GNU Scientific Library [24]: 

double F (double x, double m) 

{ return x*gsl_sf _ellint_RF(l . 0-x*x, 1.0-m*x*x, 1.0, GSL_PREC_DOUBLE) ; } 
void SNCNDN (double x, double m, double *s, double *c, double *d) 
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{ (void)gsl_sf _elljac_e(x, m, s, c, d) ; } 

• Using Numerical Recipes [25]: 

double F (double x, double m) 

{ return x*rf(1.0-x*x, 1.0-m*x*x, 1.0); } 

void SNCNDN (double x, double m, double *s , double *c, double *d) 
{ sncndn(x, 1.0-m, s, c, d) ; } 

To truly work in double precision, compared to Ref. 25, each float should be replaced by 
double and the following constants in the source code would have to be changed: 

#define ERRTOL 0.0025 
#define CA le-8 

in rf and sncndn, respectively. 
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